Install package casnet to use the functions fd_sda(), fd_psd() and fd_dfa() and plotFD_loglog(). Check the manual pages for these functions.
# Install casnet if necessary: https://github.com/FredHasselman/casnet
# !!Warning!! Beta version...
# library(devtools)
# devtools::install_github("FredHasselman/invctr")
# devtools::install_github("FredHasselman/casnet")
library(invctr)
library(casnet)
library(rio)
library(plyr)
library(dplyr)
library(ggplot2)We’ll use the data from last week:
library(rio)
library(pracma)
# Reload the data
series <- rio::import("https://github.com/FredHasselman/The-Complex-Systems-Approach-Book/raw/master/assignments/assignment_data/BasicTSA_arma/series.sav", setclass = "tbl_df")
TS1 <- rio::import("https://github.com/FredHasselman/The-Complex-Systems-Approach-Book/raw/master/assignments/assignment_data/RelativeRoughness/TS1.xlsx", col_names=FALSE)
TS2 <- rio::import("https://github.com/FredHasselman/The-Complex-Systems-Approach-Book/raw/master/assignments/assignment_data/RelativeRoughness/TS2.xlsx", col_names=FALSE)
TS3 <- rio::import("https://github.com/FredHasselman/The-Complex-Systems-Approach-Book/raw/master/assignments/assignment_data/RelativeRoughness/TS3.xlsx", col_names=FALSE)
# The Excel files did not have any column names, so let's create them in the `data.frame`
colnames(TS1) <- "TS1"
colnames(TS2) <- "TS2"
colnames(TS3) <- "TS3"
# randperm()
TS1Random <- TS1$TS1[randperm(length(TS1$TS1))]
# sample()
TS1Random <- sample(TS1$TS1, length(TS1$TS1))
TS2Random <- sample(TS2$TS2, length(TS2$TS2))
TS3Random <- sample(TS3$TS3, length(TS3$TS3))
TS3Norm <- scale(TS3$TS3)SDA:
sd based on N, not N-1, e.g. by using ts_standardise(), or use function arguments)2. If you use N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.sda on the 3 HRV time series, or any of the other series you already analysed.fd_sda to the other techniques (fd_RR,SampEn, acf, pacf)
# ACF assignment data `series`
fdS1 <- fd_sda(series$TS_1, returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 10)
## Slope = -0.64 | FD = 1.64
##
## Fit range (n = 9)
## Slope = -0.61 | FD = 1.61
##
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_sda(series$TS_2, returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 10)
## Slope = -0.58 | FD = 1.58
##
## Fit range (n = 9)
## Slope = -0.55 | FD = 1.55
##
## ~~~o~~o~~casnet~~o~~o~~~
fdS3 <- fd_sda(series$TS_3, returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 10)
## Slope = -0.64 | FD = 1.64
##
## Fit range (n = 9)
## Slope = -0.58 | FD = 1.58
##
## ~~~o~~o~~casnet~~o~~o~~~
# If you asked to return the powerlaw (returnPLAW = TRUE) the function plotFD_loglog() can make a plot
plotFD_loglog(fdS1,title = "SDA", subtitle = "Series 1")plotFD_loglog(fdS2,title = "SDA", subtitle = "Series 2")plotFD_loglog(fdS3,title = "SDA", subtitle = "Series 3")Values of other time series
# RR assignment data TS1, TS2, TS3 and shuffled versions
# This is just a check so you can see which bin sizes to expect
# It is also possible to force calclulation of specific bin sizes, see the manual pages for more info ?fd_sda
L1 <- floor(log2(length(TS1$TS1)))
L2 <- floor(log2(length(TS2$TS2)))
L3 <- floor(log2(length(TS3$TS3)))
# You can also ask for a plot directly by passing doPlot = TRUE
fdTS1 <- fd_sda(TS1$TS1[1:2^L1], doPlot = TRUE, tsName = "TS1")
##
##
## fd_sda: Sample rate was set to 1.##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.16 | FD = 1.16
##
## Fit range (n = 8)
## Slope = -0.18 | FD = 1.18
##
## ~~~o~~o~~casnet~~o~~o~~~
fdTS1r <- fd_sda(TS1Random[1:2^L1], doPlot = TRUE, tsName = "TS1 shuffled")
##
##
## fd_sda: Sample rate was set to 1.
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.67 | FD = 1.67
##
## Fit range (n = 8)
## Slope = -0.52 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
fdTS2 <- fd_sda(TS2$TS2[1:2^L2], returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.04 | FD = 1.04
##
## Fit range (n = 8)
## Slope = -0.06 | FD = 1.06
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2,subtitle = "TS2")
fdTS2r <- fd_sda(TS2Random[1:2^L2], returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.54 | FD = 1.54
##
## Fit range (n = 8)
## Slope = -0.58 | FD = 1.58
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2r,subtitle = "TS2 shuffled")
fdTS3 <- fd_sda(TS3$TS3[1:2^L3], returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.63 | FD = 1.63
##
## Fit range (n = 8)
## Slope = -0.6 | FD = 1.6
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3,subtitle = "TS3")
fdTS3r <- fd_sda(TS3Random[1:2^L3], returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.56 | FD = 1.56
##
## Fit range (n = 8)
## Slope = -0.56 | FD = 1.56
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3r,subtitle = "TS3 shuffled")
fdTS3n <- fd_sda(TS3Norm[1:2^L3], returnPLAW = TRUE)
##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.63 | FD = 1.63
##
## Fit range (n = 8)
## Slope = -0.6 | FD = 1.6
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3n,subtitle = "TS3 standardised")
# Logistic map
y1 <- growth_ac(r = 2.9, type="logistic",N = 1024)
fdL29 <- fd_sda(y1, returnPLAW = TRUE)
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 10)
## Slope = -0.41 | FD = 1.41
##
## Fit range (n = 9)
## Slope = -0.4 | FD = 1.4
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL29, subtitle = "Logistic Map: r = 2.9")
y2 <- growth_ac(r = 4, type="logistic", N = 1024)
fdL4 <- fd_sda(y2, returnPLAW = TRUE)
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 10)
## Slope = -0.64 | FD = 1.64
##
## Fit range (n = 9)
## Slope = -0.58 | FD = 1.58
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL4, subtitle = "Logistic Map: r = 4")sd based on N, not N-1)2. If you use N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.fd_psd() on the 3 HRV time series, or any of the other series you already analysed.fd_psd to the other measures (fd_RR,SampEn, acf, pacf)
# ACF assignment data `series`
fdS1 <- fd_psd(series$TS_1, doPlot = TRUE)
##
##
## fd_psd: Sample rate was set to 1.##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 512)
## Slope = 0.18 | FD = 1.57
##
## Hurvich-Deo (n = 36)
## Slope = 0.29 | FD = 1.61
##
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_psd(series$TS_2, returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 512)
## Slope = -0.27 | FD = 1.4
##
## Hurvich-Deo (n = 79)
## Slope = -0.98 | FD = 1.2
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS2,subtitle = "Series 2")
fdS3 <- fd_psd(series$TS_3, returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 512)
## Slope = 0.12 | FD = 1.54
##
## Hurvich-Deo (n = 32)
## Slope = 0.02 | FD = 1.51
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS3,subtitle = "Series 3")Values of other time series
# This is just a check so you can see which frequency range to expect
# It is also possible to force calclulation of specific bin sizes, see the manual pages for more info ?fd_psd
L1 <- floor(log2(length(TS1$TS1)))
L2 <- floor(log2(length(TS2$TS2)))
L3 <- floor(log2(length(TS3$TS3)))
# RR assignment data TS1, TS2, TS3 and shuffled versions
fdTS1 <- fd_psd(TS1$TS1[1:2^L1], doPlot = TRUE, tsName = "TS1")
##
##
## fd_psd: Sample rate was set to 1.##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.92 | FD = 1.22
##
## Hurvich-Deo (n = 33)
## Slope = -0.53 | FD = 1.32
##
## ~~~o~~o~~casnet~~o~~o~~~
fdTS1r <- fd_psd(TS1Random[1:2^L1], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.14 | FD = 1.45
##
## Hurvich-Deo (n = 21)
## Slope = -0.43 | FD = 1.35
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS1r,subtitle = "TS1 shuffled")
fdTS2 <- fd_psd(TS2$TS2[1:2^L2], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.79 | FD = 1.24
##
## Hurvich-Deo (n = 17)
## Slope = -1.12 | FD = 1.18
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2,subtitle = "TS2")
fdTS2r <- fd_psd(TS2Random[1:2^L2], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.05 | FD = 1.48
##
## Hurvich-Deo (n = 20)
## Slope = -0.01 | FD = 1.5
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2r,subtitle = "TS2 shuffled")
fdTS3 <- fd_psd(TS3$TS3[1:2^L3], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = 0.19 | FD = 1.57
##
## Hurvich-Deo (n = 16)
## Slope = 0.97 | FD = 1.79
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3,subtitle = "TS3")
fdTS3r <- fd_psd(TS3Random[1:2^L3], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.01 | FD = 1.5
##
## Hurvich-Deo (n = 27)
## Slope = 0.26 | FD = 1.6
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3r,subtitle = "TS3 shuffled")
fdTS3n <- fd_psd(TS3Norm[1:2^L3], returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = 0.19 | FD = 1.57
##
## Hurvich-Deo (n = 16)
## Slope = 0.97 | FD = 1.79
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3n,subtitle = "TS3 standardised")
# Logistic map
y1 <- growth_ac(r = 2.9, type="logistic",N = 1024)
fdL29 <- fd_psd(as.numeric(y1), returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 512)
## Slope = -1.26 | FD = 1.16
##
## Hurvich-Deo (n = 22)
## Slope = -5.23 | FD = 1.08
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL29, subtitle = "Logistic Map: r = 2.9")
y2 <- growth_ac(r = 4, type="logistic", N = 1024)
fdL4 <- fd_psd(as.numeric(y2), returnPLAW = TRUE)
##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 512)
## Slope = 0.12 | FD = 1.54
##
## Hurvich-Deo (n = 32)
## Slope = 0.02 | FD = 1.51
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL4, subtitle = "Logistic Map: r = 4")DFA analyses are not really necessary, except perhaps:
2. If you use N <- log2(length(y)), the number you get is 2^N. You need an integer power for the length in order to create equal bin sizes. You can pad the series with 0s, or make it shorter before analysis.fd_dfa() on the 3 HRV time series, or any of the other series you already analysed.fd_dfa to the other measures (fd_RR,SampEn, acf, pacf)
# ACF assignment data `series`
fdS1 <- fd_dfa(series$TS_1, doPlot = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.44 | FD = 1.55
##
## Exclude large bin sizes (n = 25)
## Slope = 0.45 | FD = 1.55
##
## ~~~o~~o~~casnet~~o~~o~~~
fdS2 <- fd_dfa(series$TS_2, returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 1.15 | FD = 1.15
##
## Exclude large bin sizes (n = 25)
## Slope = 1.31 | FD = 1.11
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS2,subtitle = "Series 2")
fdS3 <- fd_dfa(series$TS_3, returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.46 | FD = 1.54
##
## Exclude large bin sizes (n = 25)
## Slope = 0.48 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdS3,subtitle = "Series 3")Values of other time series
# This is just a check so you can see which bin sizes to expect
# It is also possible to force calclulation of specific bin sizes, see the manual pages for more info ?fd_dfa
L1 <- floor(log2(length(TS1$TS1)))
L2 <- floor(log2(length(TS2$TS2)))
L3 <- floor(log2(length(TS3$TS3)))
# RR assignment data TS1, TS2, TS3 and shuffled versions
fdTS1 <- fd_dfa(TS1$TS1[1:2^L1], doPlot = TRUE, tsName = "TS1")
##
##
## fd_dfa: Sample rate was set to 1.##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.93 | FD = 1.23
##
## Exclude large bin sizes (n = 25)
## Slope = 0.98 | FD = 1.21
##
## ~~~o~~o~~casnet~~o~~o~~~
fdTS1r <- fd_dfa(TS1Random[1:2^L1], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.58 | FD = 1.44
##
## Exclude large bin sizes (n = 25)
## Slope = 0.54 | FD = 1.47
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS1r,subtitle = "TS1 shuffled")
fdTS2 <- fd_dfa(TS2$TS2[1:2^L2], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 1.23 | FD = 1.13
##
## Exclude large bin sizes (n = 25)
## Slope = 1.31 | FD = 1.11
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2,subtitle = "TS2")
fdTS2r <- fd_dfa(TS2Random[1:2^L2], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.47 | FD = 1.53
##
## Exclude large bin sizes (n = 25)
## Slope = 0.5 | FD = 1.5
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS2r,subtitle = "TS2 shuffled")
fdTS3 <- fd_dfa(TS3$TS3[1:2^L3], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.41 | FD = 1.58
##
## Exclude large bin sizes (n = 25)
## Slope = 0.45 | FD = 1.54
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3,subtitle = "TS3")
fdTS3r <- fd_dfa(TS3Random[1:2^L3], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.47 | FD = 1.53
##
## Exclude large bin sizes (n = 25)
## Slope = 0.48 | FD = 1.51
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3r,subtitle = "TS3 shuffled")
fdTS3n <- fd_dfa(TS3Norm[1:2^L3], returnPLAW = TRUE)
##
##
## fd_dfa: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.41 | FD = 1.58
##
## Exclude large bin sizes (n = 25)
## Slope = 0.45 | FD = 1.54
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdTS3n,subtitle = "TS3 standardised")
# Logistic map
y1 <- growth_ac(r = 2.9, type="logistic",N = 1024)
fdL29 <- fd_dfa(y1, returnPLAW = TRUE)
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.48 | FD = 1.52
##
## Exclude large bin sizes (n = 25)
## Slope = 0.62 | FD = 1.41
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL29, subtitle = "Logistic Map: r = 2.9")
y2 <- growth_ac(r = 4, type="logistic", N = 1024)
fdL4 <- fd_dfa(y2, returnPLAW = TRUE)
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Detrended FLuctuation Analysis
##
## Full range (n = 30)
## Slope = 0.46 | FD = 1.54
##
## Exclude large bin sizes (n = 25)
## Slope = 0.48 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
plotFD_loglog(fdL4, subtitle = "Logistic Map: r = 4")Look at the explanation of surrogate analysis on the TiSEAN website (Here is direct link to the relevant sections)
We’ll look mostly at three different kinds of surrogates:
nonlinearTseries has a function FFTsurrogate and package fractal surrogate. Look them up and try to create a surrogate test for some of the time series of the previous assignments.
fractal::surrogate(), choose either aaft (amplitude adjusted Fourier transform) or phase (phase randomisation)nonlinearTseries::FFTsurrogates will calculate phase randomised surrogates.casnet there’s a function that will calculate and display a point-probability based on a distribution of surrogate measures and one observed measure: plotSUR_histCreate a number of surrogates, calulate the slopes and plot the results.
sda, then we do the same for psd,but you can take any of the measures as well)library(fractal)
# RR assignment data `TS1,TS2,TS3`
# Now we'll trim the time series length to a power of 2 to get the cleanest results, this is not strictly necessary
L1 <- floor(log2(length(TS1$TS1)))
L2 <- floor(log2(length(TS2$TS2)))
L3 <- floor(log2(length(TS2$TS2)))
TS1_FDsda <- fd_sda(TS1$TS1[1:2^L1])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.16 | FD = 1.16
##
## Fit range (n = 8)
## Slope = -0.18 | FD = 1.18
##
## ~~~o~~o~~casnet~~o~~o~~~
TS2_FDsda <- fd_sda(TS2$TS2[1:2^L2])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.04 | FD = 1.04
##
## Fit range (n = 8)
## Slope = -0.06 | FD = 1.06
##
## ~~~o~~o~~casnet~~o~~o~~~
TS3_FDsda <- fd_sda(TS3$TS3[1:2^L3])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.63 | FD = 1.63
##
## Fit range (n = 8)
## Slope = -0.6 | FD = 1.6
##
## ~~~o~~o~~casnet~~o~~o~~~
y1<-growth_ac(r = 2.9,type="logistic",N = 2^L1)
TSlogMapr29 <- fd_sda(y1)$fitRange$FD##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.4 | FD = 1.4
##
## Fit range (n = 8)
## Slope = -0.37 | FD = 1.37
##
## ~~~o~~o~~casnet~~o~~o~~~
y2<-growth_ac(r = 4,type="logistic", N = 2^L1)
TSlogMapr4 <- fd_sda(y2)$fitRange$FD##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.97 | FD = 1.97
##
## Fit range (n = 8)
## Slope = -0.62 | FD = 1.62
##
## ~~~o~~o~~casnet~~o~~o~~~
# These were rendomised
TS1R_FDsda <- fd_sda(TS1Random[1:2^L1])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.67 | FD = 1.67
##
## Fit range (n = 8)
## Slope = -0.52 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
TS2R_FDsda <- fd_sda(TS1Random[1:2^L2])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.67 | FD = 1.67
##
## Fit range (n = 8)
## Slope = -0.52 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
TS3R_FDsda <- fd_sda(TS1Random[1:2^L3])$fitRange$FD##
##
## fd_sda: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Standardised Dispersion Analysis
##
## Full range (n = 9)
## Slope = -0.67 | FD = 1.67
##
## Fit range (n = 8)
## Slope = -0.52 | FD = 1.52
##
## ~~~o~~o~~casnet~~o~~o~~~
fractal::surrogate()
unclass() otherwise R doesn’t recognise the objectt() trnasposes the data from 39 rows to 39 columnsfd_sda().# NOW CREATE SURROGATES
library(dplyr)
# For a two-sided test at alpha = .05 we need N=39
Nsurrogates <- 39
TS1surrogates <- as.data.frame(t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=TS1$TS1[1:2^L1] ,method="aaft")))))
colnames(TS1surrogates) <- paste0("S",1:NCOL(TS1surrogates))
# Let's look at the surrogates and originals series...
# NOTE: This is a bit elaborate, because ggplot2 wants to have data in tidy (= long) format, which means all the time series data have to in 1 column!
# Make data frame to plot:
TS1surrogates$Obs <- TS1$TS1[1:2^L1]
TS1surrogates$Time <- 1:NROW(TS1surrogates)
df <- TS1surrogates %>% tidyr::gather(key=variable, value = Y, -Time)
# This will make the observed series red
eval(parse(text = c("col <- c('Obs' = 'red3',",paste0("'",colnames(TS1surrogates)%[%39,"'"," = 'black'", collapse=","),")")))
ggplot(df,aes(x=Time, y = Y, group = variable, colour = variable)) +
geom_path(show.legend = FALSE) +
facet_grid(variable~.) +
scale_colour_manual(values = col) +
scale_x_continuous(expand = c(0,0)) +
theme_bw() +
theme(axis.text = element_blank(), strip.text = element_blank(), strip.background = element_blank())
# Now we calculate FD for each series
TS1surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TS1surrogates[,s], silent = TRUE)$fitRange$FD)
TS2surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=TS2$TS2[1:2^L2] ,method="aaft"))))
colnames(TS2surrogates) <- paste0("aaftSurr",1:NCOL(TS2surrogates))
TS2surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TS2surrogates[,s], silent = TRUE)$fitRange$FD)
TS3surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=TS3$TS3[1:2^L3] ,method="aaft"))))
colnames(TS3surrogates) <- paste0("aaftSurr",1:NCOL(TS3surrogates))
TS3surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TS3surrogates[,s], silent = TRUE)$fitRange$FD)
TSr29surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y1, method="aaft"))))
colnames(TSr29surrogates) <- paste0("aaftSurr",1:NCOL(TSr29surrogates))
TSr29surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TSr29surrogates[,s], silent = TRUE)$fitRange$FD)
TSr4surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y2, method="aaft"))))
colnames(TSr4surrogates) <- paste0("aaftSurr",1:NCOL(TSr4surrogates))
TSr4surrogates_FD <- laply(1:Nsurrogates,function(s) fd_sda(y=TSr4surrogates[,s], silent = TRUE)$fitRange$FD)x <- data.frame(Source = c("TS1", "TS2", "TS3", "TSlogMapr29", "TSlogMapr4"),
FDsda = c(TS1_FDsda, TS2_FDsda, TS3_FDsda, TSlogMapr29, TSlogMapr4),
FDsdaRandom = c(TS1R_FDsda, TS2R_FDsda, TS3R_FDsda,NA,NA),
FDsdaAAFT.median= c(median(TS1surrogates_FD),median(TS2surrogates_FD),median(TS3surrogates_FD),median(TSr29surrogates_FD),median(TSr4surrogates_FD)))plotSUR_hist() to see the results of the test.plotSUR_hist(surrogateValues = TS1surrogates_FD, observedValue = TS1_FDsda, sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS1")plotSUR_hist(surrogateValues = TS2surrogates_FD, observedValue = TS2_FDsda, sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS2")plotSUR_hist(surrogateValues = TS3surrogates_FD, observedValue = TS3_FDsda,sides = "two.sided", doPlot = TRUE, measureName = "sda FD TS3")plotSUR_hist(surrogateValues = TSr29surrogates_FD, observedValue = TSlogMapr29,sides = "two.sided", doPlot = TRUE, measureName = "sda FD logMap r=2.9")plotSUR_hist(surrogateValues = TSr4surrogates_FD, observedValue = TSlogMapr4, sides = "two.sided", doPlot = TRUE, measureName = "sda FD logMap r=4")FD: SPECTRAL SLOPE
# RR assignment data `TS1,TS2,TS3`
L1 <- floor(log2(length(TS1$TS1)))
L2 <- floor(log2(length(TS2$TS2)))
L3 <- floor(log2(length(TS2$TS2)))
TS1_FDpsd <- fd_psd(TS1$TS1[1:2^L1])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.92 | FD = 1.22
##
## Hurvich-Deo (n = 33)
## Slope = -0.53 | FD = 1.32
##
## ~~~o~~o~~casnet~~o~~o~~~
TS2_FDpsd <- fd_psd(TS2$TS2[1:2^L2])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.79 | FD = 1.24
##
## Hurvich-Deo (n = 17)
## Slope = -1.12 | FD = 1.18
##
## ~~~o~~o~~casnet~~o~~o~~~
TS3_FDpsd <- fd_psd(TS3$TS3[1:2^L3])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = 0.19 | FD = 1.57
##
## Hurvich-Deo (n = 16)
## Slope = 0.97 | FD = 1.79
##
## ~~~o~~o~~casnet~~o~~o~~~
# Logistic map - make it the length of L1
y1<-growth_ac(r = 2.9,type="logistic",N = 2^L1)
# Turn standardisation and detrending off!
TSlogMapr29 <- fd_psd(y1, standardise = FALSE, detrend = FALSE)$fitRange$FD##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.95 | FD = 1.21
##
## Hurvich-Deo (n = 13)
## Slope = -4.69 | FD = 1.08
##
## ~~~o~~o~~casnet~~o~~o~~~
y2<-growth_ac(r = 4,type="logistic", N = 2^L1)
TSlogMapr4 <- fd_psd(y2, standardise = FALSE, detrend = FALSE)$fitRange$FD##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = 0.08 | FD = 1.53
##
## Hurvich-Deo (n = 25)
## Slope = 0.44 | FD = 1.66
##
## ~~~o~~o~~casnet~~o~~o~~~
TS1R_FDpsd <- fd_psd(TS1Random[1:2^L1])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.14 | FD = 1.45
##
## Hurvich-Deo (n = 21)
## Slope = -0.43 | FD = 1.35
##
## ~~~o~~o~~casnet~~o~~o~~~
TS2R_FDpsd <- fd_psd(TS1Random[1:2^L2])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.14 | FD = 1.45
##
## Hurvich-Deo (n = 21)
## Slope = -0.43 | FD = 1.35
##
## ~~~o~~o~~casnet~~o~~o~~~
TS3R_FDpsd <- fd_psd(TS1Random[1:2^L3])$fitRange$FD##
##
## fd_psd: Sample rate was set to 1.
##
##
## ~~~o~~o~~casnet~~o~~o~~~
##
## Power Spectral Density Slope
##
## All frequencies (n = 256)
## Slope = -0.14 | FD = 1.45
##
## Hurvich-Deo (n = 21)
## Slope = -0.43 | FD = 1.35
##
## ~~~o~~o~~casnet~~o~~o~~~
fractal::surrogate()
unclass() otherwise R doesn’t recognise the objectt() trnasposes the data from 39 rows to 39 columnsfd_sda().# NOW CREATE SURROGATES
library(dplyr)
# For a two-sided test at alpha = .05 we need N=39
Nsurrogates <- 39
TS1surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(fractal::surrogate(x=TS1$TS1[1:2^L1] ,method="aaft"))))
colnames(TS1surrogates) <- paste0("aaftSurr",1:NCOL(TS1surrogates))
TS1surrogates_FD <- laply(1:Nsurrogates,function(s){fd_psd(y=TS1surrogates[,s], silent = TRUE)$fitRange$FD})
TS2surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=TS2$TS2[1:2^L2] ,method="aaft"))))
colnames(TS2surrogates) <- paste0("aaftSurr",1:NCOL(TS2surrogates))
TS2surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TS2surrogates[,s], silent = TRUE)$fitRange$FD)
TS3surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=TS3$TS3[1:2^L3] ,method="aaft"))))
colnames(TS3surrogates) <- paste0("aaftSurr",1:NCOL(TS3surrogates))
TS3surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TS3surrogates[,s], silent = TRUE)$fitRange$FD)
TSr29surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y1, method="aaft"))))
colnames(TSr29surrogates) <- paste0("aaftSurr",1:NCOL(TSr29surrogates))
TSr29surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TSr29surrogates[,s], silent = TRUE)$fitRange$FD)
TSr4surrogates <- t(ldply(1:Nsurrogates, function(s) unclass(surrogate(x=y2, method="aaft"))))
colnames(TSr4surrogates) <- paste0("aaftSurr",1:NCOL(TSr4surrogates))
TSr4surrogates_FD <- laply(1:Nsurrogates,function(s) fd_psd(y=TSr4surrogates[,s], silent = TRUE)$fitRange$FD)plotSUR_hist() to get the results of the test.plotSUR_hist(surrogateValues = TS1surrogates_FD, observedValue = TS1_FDpsd, sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS1")plotSUR_hist(surrogateValues = TS2surrogates_FD, observedValue = TS2_FDpsd,sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS2")plotSUR_hist(surrogateValues = TS3surrogates_FD, observedValue = TS3_FDpsd,sides = "two.sided", doPlot = TRUE, measureName = "psd FD TS3")plotSUR_hist(surrogateValues = TSr29surrogates_FD, observedValue = TSlogMapr29,sides = "two.sided", doPlot = TRUE, measureName = "psd FD logMap r=2.9")plotSUR_hist(surrogateValues = TSr4surrogates_FD, observedValue = TSlogMapr4, sides = "two.sided", doPlot = TRUE, measureName = "psd FD logMap r=4")x$FDpsd <- c(TS1_FDpsd,TS2_FDpsd,TS3_FDpsd,TSlogMapr29,TSlogMapr4)
x$FDpsdRandom <- c(TS1R_FDpsd, TS2R_FDpsd, TS3R_FDpsd,NA,NA)
x$FDpsdAAFT.median <- c(median(TS1surrogates_FD),median(TS2surrogates_FD),median(TS3surrogates_FD),median(TSr29surrogates_FD),median(TSr4surrogates_FD))
knitr::kable(x, digits = 2, booktabs=TRUE,formt="html")| Source | FDsda | FDsdaRandom | FDsdaAAFT.median | FDpsd | FDpsdRandom | FDpsdAAFT.median |
|---|---|---|---|---|---|---|
| TS1 | 1.18 | 1.52 | 1.24 | 1.32 | 1.35 | 1.27 |
| TS2 | 1.06 | 1.52 | 1.10 | 1.18 | 1.35 | 1.14 |
| TS3 | 1.60 | 1.52 | 1.62 | 1.79 | 1.35 | 1.58 |
| TSlogMapr29 | 1.37 | NA | 1.37 | 1.08 | NA | 1.36 |
| TSlogMapr4 | 1.62 | NA | 1.56 | 1.66 | NA | 1.56 |